package br.unb.statistic.inference;

import java.io.IOException;

import br.unb.statistic.SpectrumProcess;

public class WhittleEstimator extends HEstimator {

    private static final double tolerance = Math.pow(10, -7);
    private static final double DELTA = Math.pow(10, -9);
    private static final double SQRT_EPSILON = 1.5 * Math.pow(10, -8);
    private int numSamples;
    private double[] periodogram;

    public WhittleEstimator(Inference inference) {
        super(inference);
        this.numSamples = (int) Math.pow(2D, (int) Math.floor(Math.log(inference.getSampleSize()) / Math.log(2)));
        double[] input = new double[this.numSamples];
        try {
            inference.rewind();
            String valor = null;
            int n = 0;
            while (((valor = inference.nextValue()) != null) && (n < this.numSamples)) {
                try {
                    input[n] = Double.parseDouble(valor);
                    n++;
                } catch (NumberFormatException nfe) {
                }
            }

            inference.getDataReader().close();

            periodogram = SpectrumProcess.periodogram(input, this.numSamples);
        } catch (IOException ioe) {
        }
    }

    private double g(double H) {
        double g = 0;
        SpectrumProcess spec = new SpectrumProcess(H);
        double[] specs = spec.fgnSpectrum2(this.numSamples);
        for (int k = 1; k < this.numSamples / 2 - 1; k++) {
            g += this.periodogram[k] / specs[k];
        }
        return g;
    }

    public double getHurst() {
        double aux = (3.0d - Math.sqrt(5)) / 2;
        double beforeLastH = 0.5d + aux * 0.5;
        double beforeLastG = g(beforeLastH);
        double h = beforeLastH;
        double lastH = beforeLastH;
        double g = beforeLastG;
        double lastG = beforeLastG;
        double bottomLimit = 0.5d;
        double upperLimit = 1.0d;
        while (true) {
            double range = upperLimit - bottomLimit;
            double midpoint = (bottomLimit + upperLimit) / 2D;
            double accuracy = SQRT_EPSILON * Math.abs(h) + tolerance / 3;
            if (Math.abs(h - midpoint) + range / 2d <= 2 * accuracy) {
                return h;
            }
            double hDeviation = aux * (h >= midpoint ? bottomLimit - h : upperLimit - h);
            if (Math.abs(h - lastH) >= accuracy) {
                double aux3 = (h - lastH) * (g - beforeLastG);
                double aux4 = (h - beforeLastH) * (g - lastG);
                double aux2 = (h - beforeLastH) * aux4 - (h - lastH) * aux3;
                aux4 = 2D * (aux4 - aux3);
                if (aux4 > 0.0D) {
                    aux2 = -aux2;
                } else {
                    aux4 = -aux4;
                }
                if (Math.abs(aux2) < Math.abs(hDeviation * aux4) && aux2 > aux4 * ((bottomLimit - h) + 2D * accuracy) && aux2 < aux4 * (upperLimit - h - 2D * accuracy)) {
                    hDeviation = aux2 / aux4;
                }
            }
            if (Math.abs(hDeviation) < accuracy) {
                if (hDeviation > 0.0D) {
                    hDeviation = accuracy;
                } else {
                    hDeviation = -accuracy;
                }
            }
            double newH = h + hDeviation;
            double newG = g(newH);
            if (newG <= g) {
                if (newH < h) {
                    upperLimit = h;
                } else {
                    bottomLimit = h;
                }
                beforeLastH = lastH;
                lastH = h;
                h = newH;
                beforeLastG = lastG;
                lastG = g;
                g = newG;
            } else {
                if (newH < h) {
                    bottomLimit = newH;
                } else {
                    upperLimit = newH;
                }
                if (newG <= lastG || lastH == h) {
                    beforeLastH = lastH;
                    lastH = newH;
                    beforeLastG = lastG;
                    lastG = newG;
                } else {
                    if (newG <= beforeLastG || beforeLastH == h || beforeLastH == lastH) {
                        beforeLastH = newH;
                        beforeLastG = newG;
                    }
                }
            }
        }
    }
}
